#include <TROOT.h>

void Z_CHIRALVAL_MPS(int imom, char  *ensemble, char *pmprint)
{
//        char pmprint[8]; whichpm(ipm, pmprint);

        double ZAPREV=0.757;
        int Nmass;
        char dir[1024];

         /////////////////////////////
        // INPUT FILES///////////
        ////////////////////////////////


         // basic input file
        char input_dat[1024];
        sprintf(input_dat,"../INPUT/input_%s%s.dat", ensemble,pmprint);
      //printf("%s\n",input_dat);
        FILE *inputdat=fopen(input_dat,"r");
        fscanf(inputdat,"%i",&L);
        fscanf(inputdat,"%i",&T);
        fscanf(inputdat,"%i",&Nclust);
        fscanf(inputdat,"%s",dir);


        // mass input files

        char mu_dat[1024];
        sprintf(mu_dat,"../INPUT/mu_%s%s.dat",ensemble,pmprint);
      //printf("%s\n",mu_dat);
        FILE *mudat=fopen(mu_dat,"r");
        fscanf(mudat,"%i",&Nmass);
        double *mu=(double*)malloc(sizeof(double)*(Nmass));
        for (int im=0; im<Nmass; im++)
        fscanf(mudat,"%lf",&mu[im]);
		

       char mPCAC_dat[1024];
       sprintf(mPCAC_dat,"../INPUT/mPCAC_%s%s.dat", ensemble,pmprint);
     //printf("%s\n",mPCAC_dat);
       double mPCAC;
       FILE *mPCACdat=fopen(mPCAC_dat,"r");
       fscanf(mPCACdat,"%lf",&mPCAC);
//     printf("mpcac=%8.5f\n",mPCAC);
       fscanf(mPCACdat,"%lf",&ZAPREV);

      // total non equivalent combination of Nmass: Nmass_comb
      // I checked that Z(m1,m2)=Z(m2,m1) 
      int Nmass_comb=Nmass;
      for (int i=1; i<Nmass; i++)  Nmass_comb=Nmass_comb+(Nmass-i); 


	int aux_int;
	double aux_double;
        char MPS_dat[1024];
        sprintf(MPS_dat,"../INPUT/MPS_%s%s.dat",ensemble,pmprint);
        FILE *MPSdat=fopen(MPS_dat,"r");
        double *MPS=(double*)malloc(sizeof(double)*(Nmass_comb));
        for (int im=0; im<Nmass_comb; im++)
        fscanf(MPSdat,"%d %lf %lf",&aux_int, &MPS[im],&aux_double);

       /////////////////////////
       //// Momenta set
      //////////////////////////////
       char momentafile[1024]; sprintf(momentafile,"../INPUT/mom.dat");
      read_momenta_data(momentafile);
      Equivalent_Momenta(SinP2);



     //////////////////////////////////
     // IO files
     ///////////////////////////////

     char dir[1024];
     char file[1024];
     char z_output[1024];
     char gammaprint[8];
  
     char dir_qav[1024];

     sprintf(dir_qav,"%s/QAV/",dir);

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

  
	//////////////////////////////////
	// Read Zs
	/////////////////////////////

        double  **Zq=malloc(sizeof(double*)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust]=malloc(sizeof(double)*(Nmass));

        double  ***ZO=malloc(sizeof(double**)*(5));
        for (int igamma=0; igamma<5 ; igamma++) ZO[igamma]=malloc(sizeof(double*)*(Nclust+1));
        for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) ZO[igamma][iclust]=malloc(sizeof(double)*(Nmass_comb));

	double *aux_clust=malloc(sizeof(double)*(Nclust+1));


    // Zq

    for (int imass=0; imass<Nmass; imass++)
    {
         sprintf(file,"Zq_m%d",imass);
         read_input_onemomenta(dir_qav,file,aux_clust,imom);
	 for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust][imass]=aux_clust[iclust];
    }

    // ZO 
   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
	int imass=0;
	for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
        {
	 sprintf(file,"Z%s_m%d%d",gammaprint,im1,im2);
         read_input_onemomenta(dir_qav,file,aux_clust,imom);
         for (int iclust=0; iclust<Nclust+1; iclust++) ZO[igamma][iclust][imass]=aux_clust[iclust];
	 imass++;

       }
    }




      ///////////////////////////////////////
      // CHIRAL FIT FOR EACH MOMENTUM and PM
      //   Z= A + B*MPS
      //   where MPS^2 propto m1 + m2
      //   with m1=sqrt(mu^2+ZA^2*mpcac^2)
      // def
      //  A=ZCHI
      //  B=ZLIN
      // for ZS ZV ZT ZA and Zq
      //
      //  1/Z=A+B*MPS+C/MPS**2
      //  C=GP
      // only for ZP
     ///////////////////////////////////////


      // auxiliar for minuit routines
      double *ZAUX=malloc(sizeof(double)*(Nmass_comb));
      double **ZAUX_CLUST=malloc(sizeof(double*)*(Nmass_comb));
      for (int imass=0; imass<Nmass_comb; imass++) ZAUX_CLUST[imass]=malloc(sizeof(double)*(Nclust+1));

      double *MAUX=malloc(sizeof(double)*(Nmass_comb));

      double  **ZCHI=malloc(sizeof(double*)*(6)); // One extra for Zq
      for (int igamma=0; igamma<6 ; igamma++) ZCHI[igamma]=malloc(sizeof(double)*(Nclust+1));

      double  **invZCHI=malloc(sizeof(double*)*(6)); // One extra for Zq
      for (int igamma=0; igamma<6 ; igamma++) invZCHI[igamma]=malloc(sizeof(double)*(Nclust+1));

      double  **ZLIN=malloc(sizeof(double*)*(6)); // One extra for Zq
      for (int igamma=0; igamma<6 ; igamma++) ZLIN[igamma]=malloc(sizeof(double)*(Nclust+1));

      double  **GP=malloc(sizeof(double*)*(6)); // One extra for Zq
      for (int igamma=0; igamma<6 ; igamma++) GP[igamma]=malloc(sizeof(double)*(Nclust+1));

      double  **chi2=malloc(sizeof(double*)*(6)); // One extra for Zq
      for (int igamma=0; igamma<6 ; igamma++) chi2[igamma]=malloc(sizeof(double)*(Nclust+1));

	


    for (int iclust=0; iclust<Nclust+1; iclust++)
     {
     for (int imass=0; imass<Nmass; imass++) 
     {
                MAUX[imass]=sqrt(pow(mu[imass],2)+pow(ZAPREV*mPCAC,2));
                ZAUX[imass]=Zq[iclust][imass];
                for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=Zq[iclust_loop][imass];
     }
     lin_fit(ZAUX,ZAUX_CLUST,MAUX,Nmass,Nclust,ZCHI[5][iclust],ZLIN[5][iclust],chi2[5][iclust]); // Zq is saved in the last element (igamma=5)
     GP[5][iclust]=0.;
     sprintf(file,"Zq");  write_output_chiral_extrapolation_0GP(dir_qav,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass,ZCHI[5],ZLIN[5],chi2[5]);

     }

     // chiral fit of ZO
    // Substract Goldstone Pole for ZP
    // Simple lineal fit for the others
	

     for (int igamma=0; igamma<5; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) 
     {	
        whichGammaProm(igamma,gammaprint);
        
	if (igamma==2) // igamma=2 -> ZP
	{
		//Pseudoscalar mass square:
		for (int imass=0; imass<Nmass_comb; imass++)
		{
			MAUX[imass]=pow(MPS[imass],2);
		} 

		int imass=0;
		for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
		{
	                ZAUX[imass]=1./ZO[igamma][iclust][imass];
        	        for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=1./ZO[igamma][iclust_loop][imass];
			imass++;
		}
		lin_chiralfit_withGP(ZAUX,ZAUX_CLUST,MAUX,Nmass_comb,Nclust,invZCHI[igamma][iclust],ZLIN[igamma][iclust],GP[igamma][iclust],chi2[igamma][iclust]);
		ZCHI[igamma][iclust]=1./invZCHI[igamma][iclust];
                sprintf(file,"Z%s",gammaprint);   write_output_chiral_extrapolation_1GP(dir_qav,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass_comb,ZCHI[igamma],ZLIN[igamma],GP[igamma],chi2[igamma]);

	}
	else
	{
	        // Polar mass:
        	int imass=0;
	        for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
	        {
        	        MAUX[imass]=sqrt(pow(mu[im1],2)+pow(ZAPREV*mPCAC,2))+sqrt(pow(mu[im2],2)+pow(ZAPREV*mPCAC,2));
                	//  MAUX[imass]=mu[im1]+mu[im2];
	                imass++;
        	}

		int imass=0;
		for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
                {
	                ZAUX[imass]=ZO[igamma][iclust][imass];
        	        for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imass][iclust_loop]=ZO[igamma][iclust_loop][imass];
			imass++;
		}
		lin_fit(ZAUX,ZAUX_CLUST,MAUX,Nmass_comb,Nclust,ZCHI[igamma][iclust],ZLIN[igamma][iclust],chi2[igamma][iclust]);
		GP[igamma][iclust]=0.;
                sprintf(file,"Z%s",gammaprint);   write_output_chiral_extrapolation_0GP(dir_qav,file,imom,ZAUX,ZAUX_CLUST,MAUX,Nmass_comb,ZCHI[igamma],ZLIN[igamma],chi2[igamma]);

	}
     }

///***************************************************************************************************************///////
//********************************************************************************************************************///
///***************************************************************************************************************///////


   // Zq

   for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=ZCHI[5][iclust];
   sprintf(file,"Zq"); 
   write_output_onemomenta_real(dir_qav,file,aux_clust,SinP2_eq,imom);
  


  // ZO

   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
	for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=ZCHI[igamma][iclust];
        sprintf(file,"Z%s",gammaprint); 
        write_output_onemomenta_real(dir_qav,file,aux_clust,SinP2_eq,imom);

   } 

    delete Zq; delete ZO;
    delete aux_clust;
    delete ZAUX; delete ZAUX_CLUST; delete MAUX;
    delete ZCHI; delete invZCHI; delete ZLIN; delete GP;
  

}
